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ABSTRACT 

We use Bayesian component estimation methods to examine the prospects for large- 
scale polarized map and cosmological parameter estimation with simulated Planck 
data assuming simplified white noise properties. The sky signal is parametrized as the 
sum of the CMB, synchrotron emission, and thermal dust emission. The synchrotron 
and dust components are modelled as power-laws, with a spatially varying spectral 
index for synchrotron and a uniform index for dust. Using the Gibbs sampling tech- 
nique, we estimate the linear polarisation Q and U posterior amplitudes of the CMB, 
synchrotron and dust maps as well as the two spectral indices in ~ 4° pixels. We use 
the recovered CMB map and its covariance in an exact pixel likelihood algorithm to 
estimate the optical depth to reionization r, the tensor-to-scalar ratio r, and to con- 
struct conditional likelihood slices for Cf^ and Cf^. Given our foreground model, 
we find cr(r) ^ 0.004 for r = 0.1, cr(r) ^ 0.03 for a model with r = 0.1, and a 95% 
upper limit of r < 0.02 for r = 0.0. 



1 INTRODUCTION 

The Planck satellite mission, launched in May 2009, is mea- 
suring the polarization signal of the CMB in seven channels 
over the frequency range 30-353 GHz. It is now in normal op- 
eration and performing as expected (Planck Collaboration 
20061: rSersanelli et al.1 l2010l: iMandolesi et al.l 12010 : 
Rosset et all |201Q| : IPlanck Collaboration! 120^ . With 



a sensitivity of AT/T ~ 3 x 10~ in polarization and 
an a ngular resolution down to 5^ ("Planck Collaboration 
l2006l V Planck will produce polarization data which offer a 
multitude of opportunities including: possible recovery of 
inflationary B-modes at large scales, tighter constraints on 
the parameters describing the epoch of reionization, and 
greater understanding of the polarized nature of Galactic 
foregrounds (particularly dust, which dominates at the high 
frequency channels of Planck). 

In this paper, we adapt existing foreground separation 
tools to explore the prospects of estimating low-resolution 
polarized CMB maps from the Planck mission. We separate 
the polarization from the temperature signal and focus on 
the problem of estimating the Q and U component maps. 
From these estimated CMB maps and their covariance, we 
then use a low-^ pixel likelihood estimator to get estimates 
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on the optical depth to reionization, r, and the tensor-to- 
scalar ratio, r, and to construct low-^ conditional likelihood 
slices of Cf^ and Cf^. 

For an all-sky experiment like Planck^ component sep- 
aration of the polarization signal will be much more dif- 
ficult than for the temperature counterpart, in part, be- 
cause the ratio of the foreground signal to CMB signal 
is higher. There is a large body of work on CMB fore- 
ground su btraction for temperatu r e mea surem ents (see re- 
views by Delabrouille k, Cardosol (|2009h and I Leach et al.l 
(2008)), and a growing body of work on the problem of 
polarized foreground subtraction, including applications to 
the WMAP data and Planck simulated data. The five- 
year analysis of the WMAP data included results from 
two approaches to co mponent separation: template clean- 
ing (Gold et al. 20091) an d a Bayesian parameter estima- 
tion method ( Dunklev et al.l l2009al) whi ch we use in this 
paper. For P/ancA:. lEfetathiou et aP (|2009l ) compare an inter- 
nal linear combination technique with their template-fitting 
scheme to assess the impact of foregrounds on B-mode de- 
tection at low ^s. Betoule et al. (2009) provide constraints 
on r with Planck data and future CMB experiments us- 
ing a component sep aration method based on SMICA and 
iRicciardi et al.l (|2010l ) describe a correlated component anal- 
ysis for Planck. Finally, forecasts on polarization foreground 
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re moval for a future CM BPol mission have been presented 
in lDunklev et all (|2009bh . 

In ^2.11 we provide details of our Bayesian parameter 
technique, Gibbs sampling, and our model of the data. In 
JSl we describe the experimental specifications of Planck and 
our simulated maps. We present the estimated maps from 
the parametric technique, the estimated r and r likelihoods, 
and the low-^ conditional likelihood slices on Cf^ and Cf^ 
in SI 



2 METHODS 

In the Bayesian parameter estimation method of foreground 
separation, the emission models of the CMB and foregrounds 
are parametrized based on our understanding of their fre- 
quency dependence. We then use a sampling method to es- 
timate the marginalized CMB Q and U maps (and addition- 
ally the marginalized foreground maps) in every pixel. In this 
analysis, we use HEALPix (Gorski et al. 2005) Nside — 16 
maps containing Np = 3072 pixels. We use two different im- 
plementations of Gibbs sampling to estimate the maps and 
compare their results in the case of a simplified diagonal 
noise model. The first is a code called Commander which we 
revi ew in 2. 2 1 and refer the reader to Eriksen et al. (2006) 
and lEriksen et al.1 (|2008l ) for more details. The second is a 
code call ed Galclean wh i ch we review in ^2.31 and refer the 
reader to iDunklev et al. 1 (|2009al ) for more details. 



2.1 Bayesian Estimation of sky maps 

Consider the simple case in which the data model is given 
by 

d = s + n (1) 

where d are the observed data, consisting of Q and U polar- 
ization maps observed by Planck^ s is the sky signal, and n 
is the instrument noise. We wish to estimate the sky signal 
s which is achieved by computing the posterior distribution 
P(s|d). By Bayes' theorem, this distribution can be written 
as 

P(s|d) oc P(d|s)P(s) (2) 
with prior distribution P(s) and Gaussian likelihood 

- 21nP(d|s) = + c, (3) 

with 

x' = (d-s)*N-Hd-s) (4) 

and a normalization term c. It is straightforward to gen- 
eralize to multi-frequency data. In the case that the noise 
N is uncorrelated between channels, the likelihood can be 
written as 

x' = ;^[d,-s,]'N;i[d.-s.] (5) 

where d^^ is the observed sky map at frequency v and N^^ is 
its covariance matrix. 

We define the parametric model for the total sky signal 
in antenna temperature for our three-component model (k — 



1 for CMB, k — 2 for synchrotron emission, and k — 3 for 
thermal dust emission) as 

Sjy = ^OLk{iy;Pk)Ak (6) 

k 

where Ak are amplitude vectors of length 2Np and cxki^^] h) 
are diagonal coefficient matrices of side 2Np at each fre- 
quency. Given that the CMB radiation is blackbody and as- 
suming that the spectral index of the Galactic components 
do not vary over our frequency range, the coefficients are 
given by 

ai(!.,/?i) = aiM = /MI (7) 

Oi2{v,p2) = diB.g[(v/v3of-'] (8) 

a3(!^,/?3) =diag[(!./i^353)'^^]. (9) 

Here we have defined the function f^v) which converts the 
CMB signal I from thermodynamic to antenna temperature, 
and the two spectral index vectors (52 and /^s for synchrotron 
and dust, respectively. We also set the pivot frequencies to 
30 GHz and 353 GHz. 

Though the spectral indices for Q and U in a given pixel 
are expected to be similar (following from the assumption 
that the polarization angle does not change with frequency) , 
we allow the option for the indices to be sampled indepen- 
dently for Q and for U. Thus, our model is completely de- 
scribed by QNp amplitude parameters A = (Ai,A2,A3) 
and ANp spectral index parameters (3 — (/3^, /3^, /3|^). 
We impose a flat prior on amplitude-type parameters and 
Gaussian priors on the spectral index parameters of (52 — 
—3.0 =b 0.3 for synchrotron and /^a = 1.5 di 0.5 for dust. The 
priors we have chosen have central value and standard devi- 
ation at approximately the average and range of values typi- 
cally o bserved and sirnulated ( see, for example JPraisse et al.l 
(2008): lD^klev et al.l (|2009bl ) for further discussion). 

With our data model and priors defined, our aim is 
to estimate the joint CMB-foreground posterior P(A,/3|d) 
from which we can then obtain the marginalized distribution 
for the CMB amplitude vector as 

p(Ai,d) = J p{A,/3\d)dA2dA3d/3 (10) 

and similarly for the other model parameters. 

For the multivariate problem that we are considering, 
Gibbs sampling draws from the joint distribution by sam- 
pling each parameter conditionally as follows 

A^+' ^P(A|/3,d) (11) 

/3*+' ^P(/3|A,d) (12) 

In the next two sections, we compare and contrast the differ- 
ent methods that Commander and Calclean implement to 
sample the amplitude-type and spectral index parameters. 

2.2 Gibbs sampling with Commander 

Commander is a flexible code for joint com ponent separation 
and CMB pow er spect rum estim ation (Jewell et al. 2 0041 : 
I Wandelt et 'aD'2004: Eri ksen et al . 2004; Larson et al. 20071). 
Commander has typically been used in its full implementa- 
tion, in which the parametric model of the total sky signal is 
sampled jointly with the CMB power spectrum. We describe 
here how to use Commander to do sampling of the sky signal 
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only. The theory of Gibbs samphng allows the joint density 
P{s,C£\d) to be sampled by alternately sampling from the 
two conditional densities 



P(s|C;,d) 



Ci+^ ^ P{Ce\s*+\d). 



(13) 
(14) 



The conditional sky signal distribution can be written as 

Pis\Ce, d) oc P(d|s, Ce)Pis\Ce) (15) 

^ ^-i(d-s)'N-l(d-s)g-is*S-s (^g) 

where S and N are the signal and noise covariance matrices. 

Since we are only interested in doing only the compo- 
nent separation part, this is akin to ignoring the the P{s\Ci) 
term in the above algorithm, giving the following distribu- 
tion 



P(s|d) oc e 2 



i(d-s)*N-l(d-s) 



(17) 



For our purposes of low-^ component separation, r and r es- 
timation, and evaluation of low-^ Cf^ and Cf^ conditional 
slices, it is arguably optimal to "switch off' the d sampling 
step in Commander. This is because the d sampling step is 
time-consuming and (in the case that the CMB map has ap- 
proximately Gaussian uncertainties after marginalizing over 
foregrounds) equivalent estimates on the spectra can be ob- 
tained quickly by using an exact pixel likelihood given the 
estimated CMB map and CMB covariance. 

Alternatively, and without dropping the P{s\C£) term 
in Eq. [T5l we can say that we are conditioning on = 
for the correlated CMB component, while simultaneously 
allowing for an uncorrelated CMB component with a sep- 
arate value in each pixel. The net result is that the CMB 
amplitudes are sampled in the same way as the foreground 
amplitudes, and the d sampling step has been omitted en- 
tirely. 



2.2.1 Sampling the amplitudes 

The conditional distribution P(A|/3, d) for fixed /3 is a 6Np- 
dimensional Gaussian from which it is straightforward to 
sample. First, we define the data model as 



(18) 



The conditional distribution is 

4E[d.-E«fc(^;/3fc)Afc]^N-l[d,-E«fc(j.;/3fc)Afc] 



P(A|/3,d) oc e 



OC e 2 ^ ^ ^ ^ 



(19) 
(20) 



where A is the Wiener-filter mean given by A = Fx with 
elements 

F-i = ^ar(!.;A)N;iafe,(z.;A) (21) 
x = ^afe(z^;A)N;M.. (22) 

V 

The sampling algorithm that Commander employs 
solves the set of linear equations 



where b is the Wiener-filter mean plus random fluctuations 
(given by white noise maps w^^) 

b = ^ Oik{v\ /3fc)N-M. + ^ OLk{v\ /3fc)N-^/V.. (24) 



The solution vector A has mean A and covariance matrix 
F, and the next amplitude sample is given by A^+i = A. 

2.2.2 Sampling the spectral indices 

For fixed amplitude, the spectral index sampler in Com- 
mander is a standard inversion sampler. This algorithm first 
maps out the conditional probability distribution P(/3| A, d) 
by evaluating the likelihood (given by Eq. [5]) at 500 points 
between the upper and lower 5cr limits and then com- 
putes the corresponding cumulative probability distribution 
F{/3\A) = f^^P{/3'\A,d)d/3'. Next, a random number u 
is drawn from the uniform distribution [/[0, 1]. Thus, the 
sample from P(/3|A,d) is given by F{/3\A) = u. 

In Commander, Eg. I12lcan be iterated more than once 
in each main Gibbs loop as an inexpensive way to reduce 
correlations between consecutive samples. We allow three 
spectral index iterations for each main Gibbs iteration. 

2.3 Gibbs sampling with Calclean 

Calcl ean, the Gibbs sampling algorithm described in 
Dunkl ev et all (l2009al \ solves the same set of equations out- 
lined in ^2.11 and defines the same parametric model given 
in Eqs.[6]|9l The method that Calclean uses to do the ampli- 
tude sampling, described below in ^2.3. H is similar to that of 
Commander, with a different technique for drawing a sam- 
ple. The main difference between Calclean and Commander 
is in the spectral index sampling, where Calclean implements 
a Metropolis-Hastings algorithm as outlined in ^2.3.21 Given 
enough time to converge, both methods should give the same 
estimates. 



2.3.1 Amplitude sampling 

Calclean solves the set of linear equations given by A = Fx, 
where F and x are given by Eqs.[2T]and[22], for the Wiener- 
filtered mean A. The Calclean code then draws a Gaussian 
sample G and constructs the next amplitude sample as 



Ai+i = A + L 



(25) 



F"^A = b 



(23) 



where L is the lower cholesky decomposition of the Fisher 
matrix F — LL^. 



2.3.2 Spectral index sampling 

Calclean uses the Metr o polis- H astings algorith m (see , 
for example, Kno x et al.1 (|200lh : iLewis &: BridS (2002); 
.Dunklev et al.. (,2005 1)7 to sample the index vector. Briefly, 
a trial step Atrial is drawn from a Gaussian distribution of 
width atrial and centred on the current /3 vector. The cur- 
rent and trial vectors are used to construct the current and 
trial posteriors 

-21np(/3|A,d) = x' (26) 
where is given by Eq. \5\ which includes the full noise 
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correlation matrix. Then the ratio of the trial to current 
posterior, r, determines whether to move to the trial position 
(with probability r), or to stay at the original position (with 
probability 1-r). 

2.4 Processing the sampled distribution 

We examine the results from the Commander and Galclean 
Gibbs sampling chains by forming the mean map and covari- 
ance matrix. The mean map can be found from the marginal- 
ized distribution as 



Ak\d)AkdAk 



1 



(27) 



where we sum over the no elements in the Gibbs chain 
for the kih component of the model. We find that a typ- 
ical Gibbs chain takes roughly 100 iterations to "burn-in" 
(or converge to the equilibrium distribution) when the spec- 
tral indices are initialized at random values drawn from the 
Gaussian prior. Therefore, we remove the first 100 elements 
in each Gibbs chain before processing the samples. It should 
be noted that the marginalized distribution p(Afc|d) may 
not be an exact Gaussian, but we assume it to be to good 
enough approximation. 

Similarly, the covariance matrix for Afc, is estimated by 
summing over chain elements, as shown here for the gen- 
eral noise case, for the covariance between pixel x and y for 
component k 



Cxy,k — {^x,k^y ,k) {^x,k){^y,k) 



(28) 



= — Y,{K,k - {A^mK,k' - {Ay,u)). (29) 

For diagonal noise, the covariance in pixel x for component 
k' reduces to 

-J no 

a.fc = — - {A^,k))(A,k - (^.,fe)). (30) 

Note that while we found no = 10, 000 samples to be enough 
to ensure convergence of the chains for the case of diagonal 
noise, we expect to need an order of magnitude more sam- 
ples to construct the full covariance matrix in the case of 
correlated noise. 



2.5 Likelihood estimation 

The product of a Bayesian parametric map estimation 
method is both a CMB map and a covariance matrix (which 
can be estimated from the marginalized posterior distribu- 
tion) and together these products can be used to place con- 
straints on cosmological parameters. We compute the like- 
lihood of the estimated maps, given a theoretical angular 
power spectrum, using the method described in Page et aP 
(|2007l ). 

Summarized here, the standard likelihood is given by 



L{xn\S) ' 



exp[-|m*(5' + Ar)-^m] 



(31) 



\S + N\^/^ 

where m is the data vector containing the temperature, T, 
and polarization maps, Q and U, rip is the number of pixels 
in each map, and S and N are the signal and noise co- 
variance matrices. Under the assumption that noise in the 



temperature can be ignored - which holds at low multipoles 
where the signal dominates - the standard likelihood can be 
simplified to 



L(m|5') oc 



exp[-|m*(5'p + Arp)-^m] 



\Sp + Np\^/^ 

2 



exp[-^T*5'^^T; 



(32) 



\St\'/' 



where St is the temperature signal matrix, m = (Qp, Up) 
is the new polarization data vector given by 

^ 23 ^TE ^ 

Qp = Qp - - Tim{ + 2Yim,p -\--2Y/m,p) (33) 

^ e=2 ^£ m=-e 

23 TE ^ 

UP = UP--J2^ ^^-( + 2y€m,p+-2y;^,p) (34) 

i=2 ^ m=-i 

and 5'p is the signal matrix for the new polarization vector. 
This new form of the likelihood allows us to factorize it 
into the likelihood of temperature and polarization, with 
information about their cross-correlation preserved. 

The two cosmological parameters most closely linked 
with the large scale CMB polarization signal are the optical 
depth to reionization, r, and the tensor- to-scalar ratio, r. 
The signature of reionization is at ^ < 10 in Cf^ where the 
amplitude of the reionization signal is proportional to r^. 
The tensor-to-scalar ratio r directly scales the Cf^ power 
spectrum and is best probed at low ^'s before Cf^ due to 
lensing dominates. 

By varying only the optical depth to reionization r and 
the power spectrum amplitude (such that the temperature 
anisotropy power oX I — 220 is held constant), we can cal- 
culate the likelihood at each value of r. This allows us to 
estimate limits on the optical depth to reionization includ- 
ing foreground uncertainty. Similarly, we can vary only the 
tensor-to-scalar ratio r, fixing all other parameters at con- 
cordance values, and calculate the likelihood at each value 
of r. 

The standard likelihood estimator that we use here 
assumes that the estimated CMB map has approximately 
Gaussian uncertainties after marginalizing over foregrounds. 
In practice, marginalizing over foregrounds may induce non- 
Gaussianities and there are several ways of addressing this 
issue. One option is to modify the standard likelihood to 
include a non-Gaussian term. Alternatively, the C^'s can 
be sampled jointly with the maps in a full Bayesian frame- 
work. This type of scheme is implemented in the Comman- 
der code, in which the problem of sampling from the joint 
density P{s^Ci\d) is reduced to that of sampling from the 
two conditional densities P{s\Ci^ d) and P{Ci\s^ d). We have 
already described the sampling algorithm for the first con- 
ditional distribution P{s\Ci^d) in ^2.2.11 The second condi- 
tional distribution, P(C£|s,d), reduces to P{Ci\s) since the 
data does not further constrain the power spectrum if we 
already know the CMB sky signal. Then the distribution 
reads 



(35) 



for whic h there is a sirnple te xtbook sampling algorithm de- 
tailed in lEriksen et al.1 (|2004l ). 
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In ^4.31 we will compare foreground- marginalized Ci es- 
timates from our standard pixel-likelihood code with those 
from the Commander Gibbs sampler. 

2.6 Comparison to template cleaning 

Template cleaning is an alternative method of component 
separation that assumes that the sky at any frequency can 
be modeled as a linear sum of fixed spatial templates. In 
the regime of perfect templates and no spatial spectral in- 
dex variations, template fitting would give the optimal fore- 
ground subtraction and marginalization. We do a compar- 
ison of our results to a simple template cleaning method 
which is implemented in Commander. The data model is 
given by 

N 

d. = A + ^6,/,(z.)T, (36) 

where the first term on the right-hand side is the CMB sky 
signal, the second term represents the sum over N templates 
Tj with fixed frequency scaling fj (ly) and overall amplitude 
bj. This implementation of template cleaning assumes no 
spatial variation in the frequency scaling and is therefore 
limited in usefulness to cases where the spectral index vari- 
ations in the data are small. However, it is a fast method 
that has been successfully used for the analysis of WMAP 
data, for example. 

We use the WMAP 23 GHz map for the low- frequency 
synchrotron template and the Planck simulated 353 GHz 
map as the high-frequency dust template. This leaves six 
remaining channels of Planck (30 - 217 GHz) to be used as 
data for the template fitting. These templates are fitted to 
the data, the best fit coefficients for each component are 
found and the templates are subtracted from the map us- 
ing these coefficients in order to obtain a clean CMB map 
at each frequency. These can then be optimally combined 
through inverse variance weighting, and the likelihood com- 
puted using the method in ^2.51 This method assumes that 
the templates have no associated uncertainties; methods to 
propagate templ a te un certainty have been considered in e.g. 
lEfstathiou et al.1 (|2009l ). 



3 SIMULATED MAPS 

We use a software package called the Planck Sky Model 
(PSM, version 1.6.6) developed by the Planck Working 
Group 2 to generate our simulated CMB and synchrotron 
maps. We generate maps at the seven polarized Planck fre- 
quency channels (30, 44, 70, 100, 143, 217, and 353 GHz). 
In our analysis, we do not apply beams or smoothing to 
the data and we leave the sky unmasked in the Gibbs sam- 
pling step; these issues should be included in a more realistic 
analysis. 

We use the PSM gauss ian.cosmo option to generate 
realizations of the CMB. The PSM simulates the CMB by 
feeding a set of standard ACDM cosmological parameters 
to CAMB which produces a set of corresponding ds. A 
Gaussian random CMB temperature and polarization field 
is then drawn according to this spectra. We use PSM model 
6 for the synchrotron emission. The synchrotron Q and U 



Frequency 
GHz 


FWHM 
arcmin 


AT (Q and U) 
/iK (at nside 16) 


30 


33 


1.15 


44 


24 


1.16 


70 


14 


1.16 


100 


9.5 


0.47 


143 


7.1 


0.37 


217 


5.0 


0.61 


353 


5.0 


1.85 



Table 1. Planck specifications. The sensitivity is specified for a 
pixel of area 3.66° (area of pixel of at HEALPix nside=16, npix 
= 3072). 

emission maps are given as an extrapolation in frequency of 
the polarized 23 GHz WMAP map: 

Q.(P) = Q23(P)(J)'^^'^ (37) 

U4p) = U2s{p) (J) (38) 

The model for the spectral ind ex is taken to be model 
4 of iMiville-Deschenes et all (|2008l \ given by 

_ log(P23/^/s>$'408) . . 

log(23/0.408) ^^^^ 

where P23 is the WMAP polarization map at 23 GHz, g is 
a geometrical reduction factor (reflecting depolarization due 
to magnetic field structure), fs is the intrinsic polarization 
fraction from the co smic rays energy spe ctrum, and 5*408 is 
the 408 MHz map of lHaslam et al.l (|l982h . As an initial test, 
we generate our own simple power-law du st model, extrap- 
olating from the predicted 94 GHz map in iFinkbeiner et al.l 
(1999): S^y(p) = 5*94 (p) (^)^^ We set the dust spectral in- 
dex to Pd = 1.5 uniformly over the whole sky. Although the 
observed dust emission is known to fit better to a multi- 
component model, it is reasonable to approximate it with 
a single-component model at frequencies below 300 GHz 
where the dust polarization is likely to be dominated by 
a single component. 

We generate diagonal white noise realizations based on 
the noise levels listed in Table [J taken from the Planck 
Bluebook (jPlanck Collaboration! I2OO6I I. and scale the given 
noise levels at beam-sized pixels to the corresponding noise 
level at Uside = 16 sized pixels. Our simplified noise model 
contains no 1//- noise or other correlations, the addition of 
which would increase effective noise levels. 



4 RESULTS 

We evaluate the results from Gibbs sampling in a variety 
of ways. First, we will check that our two Gibbs sampling 
codes. Commander and Calclean, give results that agree 
both in terms of their estimated maps and also in their es- 
timates on cosmological parameters. Next, we test the low-^ 
pixel likelihood code and foreground cleaning for potential 
biases. Through a number of different test cases, we then 
look at the effect of foreground cleaning on r and r esti- 
mates. Finally, we examine the conditional likelihood slices 
of the Cf^ and Cf^ power spectra and compare the power 
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at each multipole in the Gibbs-cleaned CMB map to the 
template- cleaned CMB map. 



4.1 Comparisons and Testing 

4.1.1 Comparison between Commander and Galclean 

In this section, we compare our two Gibbs sampling codes, 
Commander and Calclean, and show that they produce com- 
parable CMB and foreground amplitude maps and spectral 
index maps for a single diagonal noise test case. Then we 
use the resulting maps from four simulations to place con- 
straints on two cosmological parameters which are relevant 
to the low-^ polarization signal: the optical depth to reion- 
ization and tensor-to-scalar ratio. We show that the esti- 
mates from Commander and Calclean on r and r agree to 
less than O.lcr in all cases. Given the robust agreement be- 
tween the two Gibbs sampling codes, and since Commander 
has the advantages of being highly parallelized and more 
flexible than Galclean^ the remainder of the results in this 
section are produced using Commander only. 

In Fig. [1] we show the input and marginalized output 
Q and U polarization maps for the CMB, the difference be- 
tween the input and output maps in standard deviations 
per pixel for the Q component, 5 = {Qin — Qout)/crQ^ and 
similarly for the U component. The marginalized amplitude 
maps are scaled to the same scale as the input CMB map. 
The marginalized output CMB maps are found to look cor- 
rect (with the addition of noise) compared to the input 
CMB maps and the normalized deviations are Gaussian dis- 
tributed with a standard deviation of one. We find that the 
Commander and Calclean marginalized output CMB maps 
agree to better than 0.2a over almost 99% of the sky. The 
small differences that we see between the marginalized out- 
put maps for the two codes can be attributed to differences 
in sampling since we find differences of the same magnitude 
between output of the same code initialized with different 
seeds. We find that the difference between the input and out- 
put is greatest in the galactic plane, where the foreground 
signal is high. However, these effects are captured in the 
marginalized errors in the maps (also shown in Fig. [T|) so the 
deviation maps do not have strong spatial dependence. The 
marginalized error maps are a useful product of parametric 
map estimation and could potentially be used to define a 
mask threshold level. Assuming these white noise proper- 
ties, we find that the signal-to-noise ratio in Planck should 
be good enough to clearly see features in the Planck CMB 
polarization maps. 

Fig. [2] shows the input and output polarization ampli- 
tude P = + U2 

maps for the synchrotron and dust 
maps, and the difference in standard deviations per pixel for 
the Q component. For the synchrotron and dust amplitude 
maps, we find that the Commander and Calclean marginal- 
ized output maps agree to better than O.lcr over almost 99% 
of the sky. 

The results from the spectral index sampling for the 
thermal dust and synchrotron components are shown in 
Fig. [3] We show the input spectral index maps (which is a 
uniform 1.5 for f5d) and the marginalized output Q index. We 
also plot the error map and difference in standard deviations 
per pixel. The bottom panels of Fig. [3] show how the effect 
of our Gaussian priors of pa = 1.5 ±0.5 and ps = —3.0 ±0.3 



varies over the sky. In regions of low signal-to-noise, away 
from the Galactic plane, the prior is stronger than the like- 
lihood and the error is driven to the prior value. In contrast, 
the likelihood dominates over the prior in regions of high 
signal-to-noise. In the dust error map, we find the error to 
be as small as 0.0026 in the Galactic plane where the like- 
lihood dominates by a factor of nearly 200. Thus, for data 
containing a spatially varying /^d, we would expect to con- 
strain those variations in the high signal-to-noise regions. 

Next, we run our low-^ pixel likelihood code on the full 
CMB±foregrounds case, using the resulting CMB maps and 
covariance matrices from both Commander and Calclean. 
The r likelihoods for four simulations are plotted in Fig. 2] 
and we find that three of the four results are consistent with 
T = 0.1 at < 2(7 (the fourth is consistent at ^ 2.2cr). We 
find constraints on r for the r = 0.1 CMB±foreground case 
with maps estimated by Commander and Calclean. These 
likelihood curves are also shown in Fig. U As with the r 
estimates, we find that the r results are consistent with r = 
0.1 at < 2(7 for all of the simulations. The estimates on 
both r and r for Commander and Calclean show that the 
two codes agree to better than O.lcr. We will discuss the 
estimated errors on r and r in ^4.21 

4.1.2 Testing the likelihood code 

To check for potential biases, we have tested the low-^ pixel 
likelihood code on foreground- free simulations. We obtain 
these simulations by co-adding the 100, 143, and 217 GHz 
simulated maps, M^, with inverse- noise weightings 

M = Y.^a- (40) 

i * 

where the co-added error, cr^ is given by 



The likelihood curves for the r = 0.1 foreground- 
free case are shown in Fig. [S] We calculate the likelihoods 
for 10 simulations and take the sum of the log-likelihoods 

( ^ —2lnLi) of the 10 simulations to obtain an average 

distribution, shown in the right-hand plot of Fig. [5l We find 
that the foreground-free case is consistent with r = 0.1 at 
1.4cr. We repeat the foreground-free analysis for r = 0.1 
and 10 simulations, as shown in Fig. [6l and take the sum 
of the log- likelihoods to obtain an average likelihood. The 
average likelihood is consistent with r = 0.1 at 1.3cr. A thor- 
ough test for biases should combine results from thousands 
of simulations, however, for our small set of simulations, the 
deviations that we see from r = 0.1 and r = 0.1 are not 
statistically significant. 

4.2 Effect of foregrounds on parameter estimation 

The likelihoods for four simulations of the r = 0.0 case 
are shown in Fig. [71 for the two cases, foreground- free 
and CMB±foregrounds. The foreground-free simulations all 
show likelihoods which are sharply peaked at r — 0.0 while 
the CMB±foreground simulations show widening of the like- 
lihoods commensurate with the increased uncertainty due to 
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Figure 1. First row: input Q CMB map (left column), Commander posterior mean output Q map (middle column), and Galclean 
posterior mean output Q map (right column). Second row: marginalized error. Commander difference in standard deviations per pixel 
(middle column), and Galclean difference in standard deviations per pixel (right column) for the Q component. Third row: input U 
CMB map (left column). Commander posterior mean output U map (middle column), and Galclean posterior mean output U map (right 
column). Fourth row: As in second row but for the U component. 



the presence of the foregrounds. This effect is summarized 
m Table E] which gives the average upper 95% cut-off limits 
on estimates of r for r = 0.0 and the average estimates on 
cr(r = 0.1) and a{T = 0.1) with and without foregrounds. 
We a lso apply the standard WMAP P06 mask (Pag e et al] 
I2OO7I I, which masks about 26% of the sky, and calculate the 
likelihood distributions for the masked case. 

For our chosen simulations and modeling, we find mini- 
mal error inflation in (Jr and ar- (Tt remains nearly constant 
at ^ 0.005 in the absence or presence of foregrounds. Gr 
increases from ^ 0.02 to ^ 0.03 with the addition of fore- 
grounds. Our limits on r = show that it is more sensitive to 
the presence of foregrounds than an estimation of an r = 0.1 



signal. We find (Tr/r = 0.32 for r = 0.1 and cFr/r = 0.05 
for T = 0.1 using Com mander. Using a Fisher matrix ap- 
proach, [Betoule^^et^l] ([2009) find values of fXr/r similar to 
ours: (Tr/r — 0.34 with foregrounds and ar/r — 0.25 with 
noise only. In anothe r Fisher matrix forecast for Planck^ 
iBaumann et al.1 (|2009l ) finds (Tr = 0.011 for r = 0.01 without 
foregrounds. 

Our pixel likelihood code can be used not only to con- 
strain parameters, but also to find the power at each mul- 
tipole in the polarized power spectra. At each multipole, 
we compute the conditional likelihood as a function of Cf^ 
and Cf^ for ^ = 2 — 7 with all other multipoles held fixed 
at the fiducial ACDM values, using the method described 
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Figure 2. Maps of the polarization amphtude P = \JQ^ for the synchrotron at 30 GHz (first row) and dust at 353 GHz (third 

row). The difference in standard deviations per pixel for the Q component (second and fourth rows) indicate that the synchrotron and 
dust maps have been recovered to the expected statistical result. 



Table 2. Average upper 95% cut-off limits on estimates of r = 0, and average estimates on cr(r = 0.1) and g( t — 0.1). We find a 
foreground-free error on r that matches the size of errors found in analogous Fisher matrix forecasts for Planck d Betoule et al.l [20091 : 
iBaumann et al.|[2009l) . The effect of foregrounds is seen to infiate the error bar in the case of r but not r. The error on r for the r = 0.1 
model is amplified by a factor of ~ 1.4, and the 95% cut-off limit on r for the r = 0.0 model is amplified by a factor of ~ 3, when 
foregrounds are included. 
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Figure 3. Spectral index, /S, input map (top), output map (second row), error map (third row), and deviation map (bottom row) for 
dust (left) and synchrotron (right). Note that in areas of low signal-to-noise, the error is driven to the prior value of 0.5 for dust and 0.3 
for synchrotron in areas of low signal-to-noise. 
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Figure 4. Likelihood distributions for r (left plot) and r (right plot) for four simulations of CMB+foregrounds with r = 0.1 and r = 0.1. 
The four different simulations are represented by the black, blue, green, and red curves. Results from Commander are shown with a solid 
line and results from Galclean are shown with a dashed line. Note that the Galclean curve for the red simulation is completely overlaid 
by the Commander curve for the red simulation (left plot), and that the Galclean curve for the green simulation is completely overlaid 
by the Commander curve for the green simulation (right plot). We find cr(r) ~ 0.004 and a(r) ~ 0.03. 
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Figure 5. r = 0.1 foreground- free case for 10 simulations. The left-hand plot shows the likelihood distributions for each of the 10 
simulations, while the plot on the right-hand side is the sum of the log-likelihoods of the 10 distributions. 
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Figure 6. r = 0.1 foreground- free case for 10 simulations. The left-hand plot shows the likelihood distributions for each of the 10 
simulations, while the plot on the right-hand side is the sum of the log-likelihoods for the 10 simuations. 
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Figure 7. The solid curves show the hkehhood distributions for r for the r = 0.0 foreground- free simulations while the dashed curves 
show the likelihood distributions from maps estimated by Commander. The widening of the likelihood curves from solid to dashed line 
is commensurate with the increased uncertainty in our estimate of r due to the presence of foregrounds. 



in iNolta et alj (|2009l l . For example, the conditional likeli- 
hood of Cf^ is f{x) oc L(d|...,C3^^,Cf^ = x,C^^...). We 
compare the power at each multipole in the Gibbs CMB 
map to the template-cleaned case (described in ^2.6p and 
to the foreground-free case, shown in Fig. [8] and Fig. [9] 
for the r = 0.1 simulation. For Cf^, the results are 
consistent between the three cases. For the Cf^ spectra, 
we find that the template-cleaned conditional slices agree 
with the foreground-free curves as well, or better than, the 
Gibbs slices, indicating that the more economical template- 
cleaning method is an effective (and fast) option for fore- 
ground removal in the case of low spectral index variation 
established in our data model. However, we argue that Gibbs 
sampling should be used instead of, or in addtition to, tem- 
plate cleaning, in order to benefit from the Gibbs feature 
that the inclusion of foreground uncertainties in the covari- 
ance matrix can be propagated to the limits on r. This effect 
appears as the infiation in the Gibbs Cf^ distributions over 
the template distributions for the r = 0.0 simulation, shown, 
in particular for = 2, 4 and 5 in Fig. [101 

In Fig. [TTl we plot the results from the Commander 
template fitting. The data points are the best-fit template 
coefficients for the dust and synchrotron emission at 30, 44, 
70, 100, 143, and 217 GHz . The dashed curves show the 
emission, in antenna units, of the thermal dust for pd — 
1.5 and of the synchrotron for /3s = —3. The curves are 
normalized to the r.m.s. values of the 23 GHz and 353 GHz 
template maps for synchrotron and dust, respectively. 

4.3 Comparison between Ci estimates from Gibbs 
sampling and pixel likelihood code 

In ^2.5l we discussed a potential issue with our standard pixel 
likelihood code in the case that the marginalized distribu- 
tions p(A|d) contain non-Gaussianities. We investigate our 
CMB marginal posteriors and do find a small level of non- 
Gaussianity particularly in regions where the foreground sig- 
nal is large. We proposed several options for addressing this 
issue in ^2.51 and in this section we show a comparison be- 
tween our standard pixel-likelihood and Gibbs sampled d 
estimates in order to assess the level of non-Gaussianity seen 
in the CMB marginal posteriors. 



We run the pixel-likelihood code to compute the condi- 
tional likelihood as a function of Cf^ and Cf^ for ^ = 2 — 6 
with all other multipoles held fixed at the fiducial ACDM 
values. We additionally marginalize over Cj^ and Cf^ 
when computing the Cf^ likelihood in order to account 
for correlations between the TT, TE, and EE components. 
We neglect correlations between i values. We run the Gibbs 
sampler ( Commander) in the mode in which the CMB power 
spectra is sampled simultaneously with the foreground com- 
ponents, as described in ^2.51 In Fig. [12] we show slices 
through the Ci distribution obtained from the Gibbs estima- 
tor compared to the pixel likelihood. We find the estimates 
from the two methods to be equivalent up to small differ- 
ences. The small discrepancies between the Gibbs and pixel 
likelihood estimates are due to the pixel likelihood code us- 
ing Nside = 8, compared to the higher resolution Nside = 16 
used for the Gibbs code. Another source of differences may 
be from i — correlations present in the Gibbs samples but 
not in the pixel likelihood which estimates slices of Ci for all 
other multipoles fixed. These results indicate that it is rea- 
sonable to approximate the foreground-marginalized CMB 
pixel amplitudes as Gaussian in the pixel-based likelihood. 



5 CONCLUSIONS 

We have investigated two independent Gibbs sampling codes 
for polarized CMB foreground separation in the case of diag- 
onal noise, and power law dust and synchrotron models. We 
have constructed the large-scale posterior CMB and fore- 
ground amplitude maps as well as the dust and synchrotron 
spectral index maps using the Planck sky model and noise 
levels, and without masking the Galactic plane. We explored 
constraints on r and r for our Planck model and found that 
our Bayesian algorithms produced results consistent with 
T = 0.1 and r = 0.1 at 1- and 2a. We find a{r = 0.1) ^ 0.004 
and cr(r = 0.1) ^ 0.03. We find a 95% cut-off limit on an 
r — detection at r < 0.02. While our specific predictions on 
(Tr and (Jr are limited by our simplified noise and data mod- 
els, we have shown that the Gibbs-estimated CMB maps and 
errors capture the additional uncertainty due to the presence 
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and removal of foregrounds, which is then translated into an 
error inflation on r and r. 

There are many interesting extensions to this work that 
can be further explored using the tools described in this pa- 
per. More realistic data models can be considered, includ- 
ing two-component dust and a synchrotron curvature model. 
Other modelling effects can be considered, such as polarized 
free- free emission and polarized spinning dust. Mismatches 
between data and model should be investigated in terms of 
the amplification factor or biases that they impart to esti- 
mates on T and r. The noise model can be extended to in- 
clude the full Q/U noise covariances and also 1// noise. The 
tools adapted and developed in this work can also be used 
to estimate parameters for small-scale versus large-scale ex- 
periments (e. g. COrE) 
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Figure 9. Conditional likelihood slices for CMB for the r = 0.1 simulation, estimated from the polarization CMB maps cleaned 

using Gibbs sampling (red), compared to the template cleaned maps (blue) and to the foreground- free case (black). The dashed vertical 
line represents the true value of Cf^. 
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Figure 10. Conditional likelihood slices for CMB for the r = 0.0 simulation, estimated from the polarization CMB maps cleaned 

using Gibbs sampling (black), and compared to the template cleaned maps (blue). We expect that Gibbs sampling is needed over template 
cleaning to provide appropriate r limits and find that Gibbs sampling does tend to inflate the distributions, particularly for i = 2,4 and 
5, in this case. 
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Figure 11. R.m.s. fluctuations spectrum (antenna temperature units) of the polarized dust and synchrotron components. The symbols 
are the best-flt template coeflicients for dust (red triangles) and synchrotron (black stars). The dashed curves represent dust for = 1-5 
(red dashed line) and synchrotron for /3s = —3 (black dashed line). The CMB fluctuations (blue line) are normalized to the r.m.s. value, 
0.24 i^K (thermodynamic), of the simulated Q and U components of the CMB map. 
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Figure 12. Estimates of Cf^ (left plot) and (right plot) computed with two different techniques. At each i value, we plot the 

maximum likelihood value (tic mark), the region where the likelihood is greater than 50% of the peak value (thick line) and the region 
where the likelihood is greater than 95% of the peak value (thin line). The black lines (right side of each pair) are estimated with 
a pixel-based likelihood code with Ngide = 8- The blue lines (left side of each pair) are estimated by Gibbs sampling the maps and 
C^s simultaneously at Ngide = 16. Note that we do not show results for C^-^ as the comparison between conditional and marginal 
distributions for ^ = 2 is not feasible using this method. 



